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Recent theoretical studies have shown that demographic stochasticity can greatly increase the 
tendency of asexually reproducing phenotypically diverse organisms to spontaneously evolve into 
localised clusters, suggesting a simple mechanism for sympatric speciation. Here we study the role 
of demographic stochasticity in a model of competing organisms subject to assortative mating. We 
find that in models with sexual reproduction, noise can also lead to the formation of phenotypic 
clusters in parameter ranges where deterministic models would lead to a homogeneous distribution. 
In some cases, noise can have a sizeable effect, rendering the deterministic modelling insufficient to 
understand the phenotypic distribution. 
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I. INTRODUCTION 

Establishing the determinants of biological diversity is 
a fundamental question in biology. A particular aspect of 
this question that has attracted a great deal of attention 
is the distribution (in genotype or phenotype space) of a 
population of interacting organisms, and to what extent 
and under what conditions clusters of similar individuals 
tend to arise CHI- A better understanding of this ques¬ 
tion could shed some light onto the process of sympatric 
speciation, whereby a ‘mother’ species splits into two or 
more other species without geographic isolation Si- 

There is a history of mathematical models proposed 
to elucidate the mechanisms whereby species (as clusters 
of phenotypically similar organisms) tend to form spon¬ 
taneously due to competition [a-[l3- Models initially 
considered asexual populations, to facilitate analytical 
tractability. It was shown that spontaneous clustering 
can occur 0, i ini, but the phenomena is somewhat 
sensitive to particular assumptions about the functional 
form of the interaction kernels used [H, [13 1 question¬ 
ing the biological relevance of the findings. Models in 
which reproduction is sexual tend to lead to spontaneous 
clustering with less restrictive assumptions 0, [l3|, but 
the biological realism of the conditions required has also 
been questioned 0] . Recently, it was shown that stochas¬ 
tic effects can greatly increase the range of parameters 
for which species are formed in asexual models Bill, 
presenting stochastic pattern formation as a novel mech¬ 
anism for speciation, which has particular biological im¬ 
plications [13 • In this paper we will study the effects of 
stochasticity in a speciation model with sexual reproduc¬ 
tion. We will show that demographic stochasticity can 
increase the parameter range in which species cluster¬ 
ing is observed and that, in some cases, noise can have 
a sizeable effect, rendering the deterministic modelling 
insufficient to understand the phenotype distribution. 

The rest of the paper is organised as follows. We will 
present the model as well as its mathematical formula¬ 
tion in Sec. im In Sec. imi we will present the analysis in 
the deterministic (large population) limit, summarising 


some previously known results and presenting some new 
ones. We will then perform the analysis of the stochastic 
effects in Sec. El highlighting the cases in which noise 
effects are largest. We will conclude with a summary 
and conclusions. Some technical details are left for two 
appendixes. 

II. DESCRIPTION OF THE MODEL 

The model consists of a population of individuals which 
reproduce sexually and die through competition with 
each other. Each individual is described by its pheno¬ 
type, which determines the extent to which the organ¬ 
ism competes with other organisms, and the likelihood 
with which it can mate with a given other organism. For 
simplicity, we will assume that the phenotype is well de¬ 
scribed by a single scalar variable which can take on all 
possible real values. A simple example in which a single 
scalar variable is the main determinant of the strength 
of competition between organisms could be the beak size 
in birds or the jaw size in lizards [T^ (which determines 
the extent to which they feed on the same resources), or 
in general, body size. The model is stochastic in that 
deaths and births are random events (which, however, 
take place with probabilities determined by the state of 
the system). There are two basic processes: 

(i) Death. The death rate (probability per unit 

time) of individual i, di, is given by di = 
[Kf{xi)]~^ ~ where g{x) is the com¬ 

petition kernel that quantifies the strength with 
which two individuals with phenotype distance x 
compete, AT is a constant that controls the overall 
carrying capacity of the ecosystem and f{x) is a 
function that determines the relative intrinsic ad¬ 
vantage of phenotype x. 

(ii) Reproduction. Each individual reproduces at a rate 
one. The probability that individual i mates with 
individual j is proportional to m{xi — Xj), with m 
a function determining the strength of assortment 
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(preference for individuals which are alike). If indi¬ 
viduals i and j mate, an offspring is generated with 
phenotype given by Xq = {xi -\- Xj)/2 + C, where C, 
is a random variable with probability density func¬ 
tion r{x). This models mutation about the average 
of the parents phenotype {xi + Xj)l2. 


Xj)lK f (j){y)m(xi — y)dy = m(xi — Xj)/K(j)*m(xi). The 
probability density that their offspring has phenotype x 
is given by r(x— {xi +Xj)/2), where r{x) is the Gaussian 
probability density with variance cr^. With this in mind, 
we see that the rate density at which a new individual is 
created with phenotype x, /3{x, (p), is: 


This simple reproduction rule is justified when the 
character under consideration is determined by a large 
number of additive genes (that is, without dominance or 
epistasis [l^); in this limit, the ‘reproduction noise’, C,, 
is a Gaussian random variable with zero mean and vari¬ 
ance given by that, for simplicity, will be assumed to 
be phenotype-independent (the influence of this quantity 
will be one of the main aspects of our investigation). The 
function f{x) modulates the relative advantage of the 
phenotypes; it typically decays to zero for large |a:|, con¬ 
straining the phenotypes to a given region of interest. We 
are also assuming that the organisms are hermaphroditic, 
but this assumption will be relaxed later. The model is 
a generalisation of the asexual model studied in (TB |. 

Following [l^ , we will describe the state of the system 
by the density in phenotype space, K(j){x), where 

, N(t) 

4>{.x) = — '^5{x-Xi). ( 1 ) 

i=l 

Here Nit) is the number of organisms at time t and (5(a:) 
is the Dirac delta function. We have introduced the car¬ 
rying capacity, K, to obtain a function (j){x) that has a 
well-defined K ^ oo limit. 

When an organism with phenotype y dies, 4>{x) is mod¬ 
ified by the subtraction of delta function centred at y. 
Similarly, if a new organism with phenotype y is born, 
(j){x) is modified by the addition of a delta function at y. 
With this in mind, we define the operators by their 
action on a generic functional F[(j){x)] as: 

A^F[Pix)]=F[Pix)±^6ix-y)]. ( 2 ) 

Now suppose that j{x, p) is the density rate at which an 
individual with phenotype x dies, so that ^{x, (j))dxdt is 
the probability that an individual with phenotype in the 
interval (x,x + dx) dies in the time interval (t,t + dt), 
given that the state of the system is given by (p at time t. 
The definition of the process implies that 7 ( 0 , t) is given 
by: 


(p) 


’^W) / 

K-^(p* g(,x). 

fix) 


X — Xi) 

y)dy 


( 3 ) 


The probability that individual i, if it mates, does so 
with individual j is m{xi — Xj)/ ixiixi — Xk) = m{xi — 


N{t) _ 

Pix, (p)=J 2 ( 4 ) 

i,t=l 

(piz)m{y - z) 


= K j (piy) J ■ 


(p * m{y) 


K(p * m{xi) 
'-r{x - iy + z)/2)dydz. 


Combining the two contributions to the change in (p, the 
probability density of finding the system at state (p at 
time t, Pip, t), changes in time according to the following 
functional master equation [isj ]: 


^PiP,t) = J[iA^ -l)/3{p,x)P{p,t) 

+ i^t - ^hiP,x)P{p,t)]dx. (5) 


Due to the non-linearity of the system (that arises due 
to the interactions) we are unable to obtain an exact so¬ 
lution and some approximation scheme is needed to pro¬ 
ceed. Expanding the operators in Eq. to second 
order in K~^, we can derive a functional Fokker-Planck 
equation, given in Appendix For our purposes, it is 
clearer to work with the equivalent stochastic differential 
equation which takes the form (see Eqs. (IA4I) and (IA5I1 1 


dp{x, t) 
dt 

f Piy,t)Piz,t) 
J p* m{y) 
r]ix,t) 


= - j Pix,t)piy,t) 


aix - y) 


dy 


fix) 

m{y — z)r(x — iy + z)/2) dydz 


Vk ’ 


( 6 ) 


where r]{x, t) is a Gaussian white noise with zero mean 
and with a correlator which is given by Eq. (IA5I1 of the 
Appendix. It is interesting to compare this equation with 
the analogous equation in the asexual case. There m{x — 
y) = 5{x — y) and so p * m{x) = p{x,t). Then Eq. ([6]) 
becomes 

+ J Piy, t)r{x -y)dy+ , (7) 

which is the equation found in Ref. [l5l| . apart from the 
function f{x), which was not included in the form of 
the model previously analysed (note that we are allowing 
self-fertilization since, for simplicity, in Eq. o we do 
not exclude i = j] forbidding the i = j case would add 
0{1/K) terms). We will now analyse Eq. (jS]), first of 
all in the deterministic limit, and then in the general 
stochastic setting. 
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III. DETERMINISTIC ANALYSIS 

The deterministic limit corresponds to taking A —>■ oo, 
and so the governing equation is simply Eq. (jS]), but 
with the last (noise) term absent. Some progress may 
be made analytically if we assume that the ecological 
functions, namely the competition kernel [ 5 ( 2 :)], the mat¬ 
ing function [m(a;)], the function modulating the carry¬ 
ing capacity [/(a;)] and the offspring distribution [r(a;)], 
are all Gaussian functions. We denote their variances 
by and cr^, respectively. In this case, the de¬ 

terministic equation has a Gaussian stationary solution, 
= Ce~^ with CTgj and C both satisfying 

complicated algebraic equations; the ^nation for be¬ 
ing first derived by Doebeli et. al. These authors 

also found that numerical integration of the determin¬ 
istic equation showed that when cr^ and cr^ are small 
compared with a^, there is an intermediate range of cr^, 
^ ~ ~ CTm j which the Gaussian solution 

becomes unstable and a multi-modal stationary solution 
is obtained. 

The random mating case (cr^ —>■ 00 ) is particularly 
interesting to analyse. In this limit, the equation for 
reduces to a cubic equation + a4tTst + 0!20'st + Q^o = 0 , 
where 


a4 = 2a^ + a^, a2 = - 2a^(2aj - a^), 

ao = -2a^a^o-j. ( 8 ) 

Since 04 > 0 and ag < 0, the cubic equation always 
has a single positive solution. This then is the required 
solution. In this limit, the constant C is found to be 


C = 



(9) 


(we have not taken f(x) to be normalised; f(x) = 
exp(—x^/(2aj)), so that tr^ controls the size of pheno¬ 
type space available). If we now, in addition, investigate 
the > 00 limit (the most relevant when the phenotype 
distribution is not too constrained by the external fitness 
landscape), we find two very different regimes, depending 
on the relative width of the competition kernel and the 
reproduction noise distribution. When the reproduction 
noise, cr^, is smaller than a critical value given by 0-^/4:, 
the variance of the phenotype distribution is small and in¬ 
dependent of (jj, = 2a^a^/((7^ — 4cr^). If depends 
on tTj in a way which means that it diverges as Cj —>■ 00 , 
then the term dominates over the term quartic in Cst 
and the term quadratic in CTst dominates over ag- There¬ 
fore, (Tst ~ ^a‘j{4(T'^ — cr^), and this implies that if the 

reproduction noise exceeds the critical value, cr^ > cr^/4 
then the variance of the phenotype distribution grows lin¬ 
early with Uf. In both cases, numerical integration of the 
deterministic equation shows that the Gaussian distribu¬ 
tion is always stable and the phenotype distribution is 
always uni-modal, a consequence of the random mating. 


The factors determining the transition into a multi¬ 
modal distribution can be understood in a simpler way if 
we assume that the range of phenotype space is finite, for 
instance — tt < 2 : < tt, and assume that the deterministic 
equation (Eq. ([ 6 ]) without the noise term) satisfies peri¬ 
odic boundary conditions. The competition function, g, 
the assortment function, m, and the offspring distribu¬ 
tion will be assumed to be periodic functions of period 
2ty. Since there is now no need for the function f{x) 
to regulate the competition process at large lx], we set 
/(x) = I. While the periodic boundary condition as¬ 
sumption is biologically unrealistic, we expect it to have 
a small impact when the scales of the competition, mat¬ 
ing and offspring distributions are all much smaller than 
the available region of phenotype space (as determined 
by/). 

Under these conditions, the deterministic equation has 
a uniform solution (j){x, t) = constant. Using the normal¬ 
isation g{x) dx = 1, one finds that (f>st = 1. A linear 
stability analysis of the deterministic equation around 
the solution ^ = I is performed in Appendix This 
shows that the uniform solution is unstable if: 


2rfem_fe/2 - rfcmfcm_fe/2 - 1 - gt > 0, (10) 


for some fc, with ri~,rnk and gk the Fourier modes of 
the reproduction noise, the mating and the competition 
kernel, respectively. An equivalent expression to m in 
a slightly different model was derived in Q. 

From the inequality one can obtain the bifurcation 
diagram of the system that shows the regions of parame¬ 
ter space in which the deterministic equation [Eq. (IA 6 I) ] 
leads to patterns. Figure [T] portrays several bifurcation 
diagrams for the case in which the ecological functions 
are uniform and Gaussian functions (projected onto the 
(ct^jCt^) plane). The figure shows that phenotypic clus¬ 
ters are observed for low values of the reproduction noise 
am and for intermediate values of the assortativity scale 
CTa, in accordance with what was found in [l3 |. Intu¬ 
itively, clusters tend to form with a distance, d, between 
them so that they barely compete (so d ^ CTc), but these 
clusters will only be durable if individuals from different 
clusters do not mate (so d > aa), suggesting that clus¬ 
ters will not form if Ua^ Oc- Sexual reproduction tends 
to concentrate the phenotype distribution, so a moder¬ 
ate value of a a can promote the formation of clusters. 
Moreover, the clusters can be well-defined only if the re¬ 
production noise is not too large (om < d). 

In the asexual case BO it was found that the range 
of parameters for which patterns occur was much greater 
than that predicted by a stability analysis of the kind 
carried out above. So motivated by the expectation that 
the same may be true in the case of sexual reproduction, 
we go on to analyse our model for finite values of AT, 
when stochastic fluctuations will be present. 
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FIG. 1: Region of stability of the homogeneous solution of 
the deterministic equation in space (mating range - 

mutation range), for several values of (competition range), 
for uniform (upper panel) and Gaussian (lower panel) forms 
of the ecological functions. The homogeneous solution is un¬ 
stable below the lines, leading to the appearance of patterns. 
The kinks are produced when the maximum value of the left- 
hand side of the inequality changes from one value of k 
to another. The black dot marks parameter values used in 

Fig-H 


IV. STOCHASTIC EFFECTS 
A. Weak noise effects 

Numerical simulations of the individual-based model 
show that, for moderate values of the carrying capacity, 
K, the observed phenotype distribution does not always 
agree with the results obtained from the deterministic 
equation (IA6I) . For example, at points in parameter space 
corresponding to a stable uniform solution, but not too 
far from the instability boundary, one can observe clus¬ 
ters in phenotype space forming for small values of K 
(see Fig.©. 

The origin of these stochastic patterns can be under¬ 
stood by analysing the stochastic differential equation 


FIG. 2: Time evolution of the phenotype distribution for = 
0.3, = 0.35, (T^ = 0.043 and uniform ecological functions, 

in the region of deterministic stability of the homogeneous 
solution (marked as a black dot in Fig. [TJ. The carrying 
capacity is A = 20 (upper panel) and K = 200 (lower panel). 
Note that for the smaller carrying capacity clear phenotypic 
clusters form. 

Specifically, we apply the van Kampen system- 
size expansion [20| by expanding (j> about the homoge¬ 
neous solution found in Sec. Ell and writing 4‘{x,t) = 
1 The factor reflects the nature of 

the fluctuations at large A, and ^{x,t) is a new stochas¬ 
tic field. The bulk of the calculation is exactly the same 
as that carried out when performing the linear stabil¬ 
ity analysis in the deterministic case, except in this case 
the existence of the factor ensures that the noise 

term in Eq. m is retained. Thus, going over to Fourier 
variables and using Eq. (IA8I) . one finds that 

= [27’fem_fc/2 - 1 - fffc - rfemfcm_/c/2] ^fc(<) 

+ Vk{t), (11) 

where %(!) is the Fourier transform of r](x^ t). The corre¬ 
lation function of this noise is only calculated to leading 
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order within the linear noise approximation, that is, set¬ 
ting (j){x,t) = 1. This gives, using Eqs. (IA2L (IA3I) and 
(IA5I1 . f X, y) dy = 2, and therefore 

{r]k{t)v-k{t')) = 2 (27r) S{t - t'), (12) 

showing that, as usual, the noise is additive within the 
linear noise approximation. 

The size of the stochastic patterns can be quantified 
by looking at the spatial covariance of the phenotype 
distribution in the stationary state: 

Coy{4'{x), (j){x + A)) = 

{{(l){x) - {(t){x))) {(j){x + A)- {(j){x + A)))). (13) 

If the distribution in phenotype space shows high-density 
regions separated by low-density ones with a well-defined 
average distance, the spatial covariance will display a siz¬ 
able spatial modulation in A. In the linear noise regime, 
the covariance becomes 

J Cov{(j){x),4'{x + A))dx = — J {^{x)^{x + A))dx 

= («) 

n 

Using Eq. (IA13I1 . and assuming that gn,rnn and Xn are 
all even functions of n (which follows if the ecological 
functions are symmetric), we may find (^„^_„), and then 
finally obtain: 


J Cov{(j){x), (j){x + A))dx 



n= — oo 


„inA 


[l + + »"nW„TO„/2 - 2r„TO„/2] 


.(15) 



FIG. 3: Spatial covariance of the phenotype density distribu¬ 
tion. The ecological functions are uniform distributions with 
variances given by = 0.3, = 0.35, cr^ = 0.043 (marked 

as a black dot in Fig. [T|. Numerical results (symbols) were 
averaged over at least 100 measurements, after a transient of 
t = 2000. The lines correspond to the theoretical expression, 
Eq. (fT5l). 


zone close to the deterministic transition line. We will 
now show that stochastic effects can also play a more 
prominent role in the sexual case. 


B. Strong noise effects 


Expression m is compared against numerical simula¬ 
tions in Fig. |3] for the case in which the ecological func¬ 
tions are uniform distributions. For low values of the 
carrying capacity, K, the covariance shows a clear spatial 
modulation, revealing the presence of stochastic patterns. 
The theoretical expression agrees qualitatively, and the 
quantitative agreement becomes better as K increases, 
where the linear noise assumption is expected to be a 
better approximation to the stochastic dynamics. The re¬ 
sults are very similar when the ecological functions have 
other forms, for instance if they are Gaussians or belong 
to the family exp(—|a;|V( 2 cr*)) with varying 1. 

In summary, the linear noise approximation shows that 
stochastic pattern formation can lead to the formation 
of clusters in phenotype space in situations when the de¬ 
terministic description predicts a uniform distribution. 
Interestingly, the stochastic patterns (as opposed to the 
deterministic ones) decrease with the carrying capacity 
(that controls the abundance of individuals), which leads 
to particular biological implications [13 when this mech¬ 
anism is dominant. 

When comparing with the asexual case [T^, the 
stochastic patterns now appear on a relatively narrow 


A stronger noise-induced effect occurs when the mat¬ 
ing is random. As discussed in Sec. IIIIl in the random 
mating case there are two regimes: one leading to a nar¬ 
row phenotype distribution, for low reproduction noise 
( 4 ( 7 ^ < CTg), and one leading to a broad phenotype distri¬ 
bution in the large reproduction noise case. In the broad 
phenotype regime (dcr^ > simulations show that the 
phenotype distribution observed is much narrower that 
the one predicted by the deterministic analysis, see Fig. ID 
This effect happens rather generically in this regime. 

In the random mating case, since the mating range is 
the largest scale of the system, assuming periodic bound¬ 
ary conditions gives results that are quite different from 
those of the system with open boundaries. In this last 
case there is a strong bias towards the centre of the phe¬ 
notypic space. For this reason we focus on the open 
boundary conditions case, which is the more biologically 
relevant. This, however, greatly complicates the math¬ 
ematical analysis of the noise effects. Simulations show 
(Fig. m that stochastic effects lead to the formation of 
a phenotypic cluster with a width that approaches the 
deterministic prediction only for very large values of the 
carrying capacity, K. We, therefore, find that in this 
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FIG. 4: Upper panel shows the observed instantaneous pheno¬ 
type distribution in the steady state for different values of the 
carrying capacity, together with the deterministic prediction. 
Lower panel shows the evolution of the phenotype distribu¬ 
tion for K = 400, starting with the stationary distribution 
predicted by the deterministic analysis. Parameter values are 
al = 0.02, cri = 0.04, a) = 6400. 




regime demographic stochasticity has a quite large ef¬ 
fect, again leading to the formation of tight phenotypic 
clusters when the deterministic analysis predicts a broad 
distribution. Similar results are obtained when the eco¬ 
logical functions do not have a Gaussian form, showing 
the robustness of this phenomenon. 


V. THE CASE WHEN INDIVIDUALS BELONG 
TO ONE OF TWO SEXES 


So far we have assumed, for simplicity, that any two or¬ 
ganisms can mate, i.e. the organisms are hermaphroditic. 
We now go on to model the situation with two explicitly 
different sexes. 

We will denote the phenotype of female organism by 
Xi and that of the male organism by ya ■ We use different 
indices, since the number of male and female organisms 
at a given time will typically be different, and so the 


range of these indices will be different. The model is as 
in Sec. m but with the following modifications: 

(a) Only females reproduce and at rate one. 

(b) The probability that female organism i mates with 
male organism a is proportional to m{xi — ya)- 

(c) The probability that the offspring is male or female 
is 1/2. 

(d) The competition is independent of sex, so that the 

death rate of organism i is di — ~ ^j) + 

Y.a9{Xi - ya)- 

With these assumptions, one can follow through the 
analysis given in Sec.lIIIland find the analogues of Eq. ([6]): 


^ = - J (l>f{x,t)[(l)f{y,t) + ^rniy.t)] 

1 f (t>f{y,t)(l}m{z,t) 


gjx - y) 
fix) 


dy 


J 4>m* m{y) 

Vfix,t) 

Vk ’ 


m{y — z)r{x — {y + z)/2) dydz 

(16) 


and 


J (t)mix,t) [(t)f{y,t) + ())^{y,t)] 


1 r (t)f{y,t)(j)miz,t) 

2 J 4>m* m{y) 

V-mix,t) 
y/K ’ 


m{y — z)r{x — (y + z)/2) dydz 

(17) 


where the subscripts / and m denote female and male 
respectively. It is convenient to work with the sum and 
differences of and cpm'- S{x,t) = 4>f{x,t) -|- (j)mix,t) 
and D{x,t) = (j)f{x,t) — (j)rn{x,t)- If we additionally take 
A —>■ oo to obtain the deterministic equations, we find 
that 


dS{x, t) 


dt 


= -1 S{x,t)S{y,t) 


gjx - y) 
fix) 

[Siy,t) -t D{y,t)][S{z,t) - D{z,t)] 


dy 


2 J [S — D]* m{y) 

X [m{y - z)r{x - iy + z)l2) dydz], 


(18) 


and 


dD{x,t) 


= - J D{x,t)S{y,t)^^yj^ dy- (19) 


From Eq. (USD we see that a steady state solution is 
D{x) = 0, that is, (j)fix) = (j)mix)- Then the de¬ 
terministic equations for (j)f and (pm collapse into each 
other, and agree with the deterministic equation in the 
hermaphroditic case, apart from the factor of 1/2 seen in 
Eqs. (HU) and (II3- If, as in Sec. imi we take the ecological 
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functions to be all Gaussians, we again find a Gaussian 
stationary solution (/)(x)st = Ce~^ with and 

C satisfying the same complicated algebraic equations, 
except that C takes on a value of one quarter the value 
found in Sec. IIIII Note that the total population is now 
half that obtained in the hermaphroditic case because 
now we assume that only females can initiate reproduc¬ 
tion events. 

To examine the instability leading to the appearance 
of patterns, we can again assume a finite interval — tt < 
X < TT for phenotypic space, with periodic boundary con¬ 
ditions. We can now look for homogeneous stationary 
solutions, and study the stability of these solutions. It is 
clear that once again T) = 0 is a stationary solution, with 
4'f = 4’m. = 1/4 under the same conditions discussed in 
Sec. IIIII Linearising about these homogeneous solutions, 
we write S{x,t) = 4 -|- S{x,t) and D{x,t) = D{x,t). 
Keeping only linear terms in S and D, we have from 
Eq. dm) that 

= -\D{x,t) ^ D{x,t) = D{x,0)e~*^^, 

( 20 ) 

which shows that the solution D{x,t) = 0 is always sta¬ 
ble. For S'(x,t), it is more convenient to work in Fourier 
space (see Appendix 1X1) . We find that 

dSk{t) 1 r 1 ~ \ 

= 2[^fkm_k/2-rkmkm_k/2-^-gk\Sk{t) 

+ PkDk{t), ( 21 ) 

where Fk is a function of k only. Equation (I20p shows 
that Dk{t) decays exponentially with t, so Eq. (1^ im¬ 
plies that one obtains the same stability condition as the 
one found for in the hermaphroditic case, Eq. (nni. 
We can therefore conclude that the stability boundaries 
in the case of two sexes are identical to those found in 
the hermaphroditic case. This is confirmed by numerical 
simulations of the stochastic version of the model. We 
also find that stochastic pattern formation takes place as 
in the hermaphroditic case, with results from the case of 
two sexes being equivalent to those of the hermaphroditic 
case, but with a factor of 1/4 in the carrying capacity. 

VI. CONCLUSION 

The precise definition of exactly what constitutes a 
species is still open to debate . One of several dif¬ 

ferent alternatives is the ‘phenotypic clustering species 
concept’, in which species correspond to distinct pheno¬ 
typic clusters, analogous to Mallet’s ‘genotypic clustering 
species concept’ [2l|. In this view, the clustering of indi¬ 
viduals in trait or gene space that we recognise as species 
is a pattern that emerges from underlying ecological and 
evolutionary mechanisms. Just as mixtures of chemical 
constituents which react and diffuse may create patterns 
(for instance, spots and stripes), so individuals which 


react (for example, compete) and diffuse (for example, 
mutate in trait or gene space) may create patterns (clus¬ 
ters). The traditional approach of the theoretical physi¬ 
cist would then be to construct a simple model to see if 
the effect appears, and if so, then see if a deeper under¬ 
standing of the effect can be gained from an analysis of 
the model. 

This is the approach that we have adopted here. We 
have started from a simple model that only contained 
birth, death (through competition) and mutation, and 
asked under which conditions clusters of individuals in 
phenotype space were formed. As discussed in the Intro¬ 
duction, this is a question which has been investigated by 
several authors, however our study focused on individu¬ 
als who gave birth only after mating with another in¬ 
dividual, whereas most previous investigations assumed 
asexual reproduction. 

We also investigated stochastic pattern formation [H, 
[ 2 ^ . Most previous work was carried out in the case of in¬ 
finite carrying capacity (in our notation, AT —>■ 00 ) where 
the governing equations are deterministic. The standard 
way of proceeding in this case is to look under what con¬ 
ditions the constant density (homogeneous) solution of 
this equation is unstable. This is carried out by perform¬ 
ing a linear stability analysis about this homogeneous 
solution. However it has been found that frequently pat¬ 
terns are still found in regions of parameter space where 
the homogeneous solution of the deterministic equation 
is stable. These patterns typically are stochastic, and 
can be found by analysing the governing equations at fi¬ 
nite carrying capacity, K. Interestingly, the scaling of 
these patterns with K has particular biological implica¬ 
tions fr^ . 

The search for stochastic patterns in models of asexual 
reproduction has been carried out previously [H, [l^ . It 
was found that noise originating in the discrete nature 
of individuals can lead to the spontaneous formation of 
species in situations where this would not happen de¬ 
terministically. The main purpose of this paper was to 
extend this to the case of sexual reproduction. We first 
supposed that individuals played both the male and fe¬ 
male role (i.e. were hermaphrodite). We found that when 
mating is assortative (i.e. organisms show preference for 
like individuals) stochastic patterns can appear in the re¬ 
gion of stability of the deterministic homogeneous solu¬ 
tion. These patterns were, however, somewhat restricted 
to parameter values not too far from the deterministic in¬ 
stability boundary. When mating is random the stochas¬ 
tic effects are stronger, and the phenotype distribution 
for moderate K is always relatively narrow, in contrast 
with the deterministic predictions. The case where the 
individuals were either male or female led to very similar 
results, which in some cases could be mapped directly 
onto the the hermaphroditic case. 

There are several extensions of the current work which 
could be carried out. One could distinguish between “ge¬ 
netic noise” (arising from recombination and mutations), 
which affects the inheritable traits, and “environmental” 




or “developmental” noise, which leads to two individu¬ 
als with same genotype to have different phenotypes and 
which is not inherited. These two types of noise are likely 
to have rather different effects in the phenotype distribu¬ 
tion. Also, if a dominant part of the genetic noise is due 
to recombination, that is, the trait considered is deter¬ 
mined by the effect of many genes in a diploid organism, 
and the effect of the different genes is additive (no epis- 
tasis or dominance), then the noise should depend on the 
position in genotype space (i.e. it would be multiplica¬ 
tive noise). Multiplicative noise would break the spa¬ 
tial symmetry and could lead to interesting effects, but 
would be more difficult to study analytically. Another 
possible extension is to include the Allee effect, which we 
have here ignored for simplicity. We believe, however, 
that the work discussed here shows that the inclusion of 
stochastic effects is vital if we wish to predict the range 
of parameters for which patterns, and therefore possibly 
species, occur. 
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Appendix A: The mesoscopic evolution equation 


In the limit of large carrying capacity, the master equa¬ 
tion ([S]) can be expanded in powers of K~^ to give the 
following functional Fokker-Planck equation (analogous 
to the derivation in Ref. [l^ for the asexual case): 

^P(<(),t) = -J J [A{(j), x,y)P{cl))] dxdy {Al) 

* in j 

where terms of K~^ and higher in the expansion have 
been neglected. Here 

- 4 ((/>, X, y) = X, y) - x, y), 

B{(l),x,y) = 'i>{(j),x,y) + E{(j),x,y), (A2) 

where 

'^{x,y) = [ rn{y - z)r{x - {y + z)l2) dz, 

E{x,y) = (i)[x)(l)(y f^^ (A3) 

A completely equivalent way of expressing the stochas¬ 
tic dynamics of the system, is to write down the equiv¬ 
alent stochastic differential equation. This takes the 
form 


d(l){x, t) 


= J -4(<?i', X, y) dy + (A4) 
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where 77 ( 0 ;, t) is a Gaussian white noise with zero mean 
and with correlator 

{r]{x,t)v{x',t')) = J B{^,x,y)dyS{x-x')6{t-t'), (A5) 

understood in the sense of Ito. 

As described in Sec. uni of the main text, we can ob¬ 
tain some insight into the transition into a multi-modal 
distribution by going over to a finite phenotypic space; 
specifically we assume that — tt < x < tt. We begin by 
analysing the deterministic dynamics of the model, which 
is found by letting A —>■ 00 in Eq. That is, 

=-y (A6) 

We substitute 


(p{x,t) = i (A7) 

into Eq. (IMI) and only keep linear terms in (j){x,t)- As¬ 
suming that (/,m and r are periodic with period 27r and 
normalised to unity in the interval (—tt, tt), and f{x) = 1, 
one finds that 

^^fc(t) = [-1 - gk + 2 rfem_/c/2 - rfcTO/cTO_fe/2] 4>kit). 

m 

Here we have gone over to Fourier space, since the linear 
nature of the problem, and the translational invariance, 
make this a natural choice. The Fourier modes are de¬ 
fined by 

hk= f h{x)e-^^^dx, h{x) = — y (A9) 

J- 2^ k 

From Eq. (IA8I) we see that if the condition given in 
Eq. (fTOl) of the main text holds, then the homogeneous 
solution (/) = 1 is unstable. 

If we wish to carry out a system-size expansion, as 
discussed in Sec. IIV A[ then we write (j){x,t) = 1 + 
K~^/^^{x, t) and expand in terms of As discussed 

in the main text, this leads to Eq. (HU: 

^&(0 = -Pkik{t) + r]k(t), (AlO) 

where pfc = l+gk-‘2rkm_k/2+rk'mkm_k/2- Multiplying 
by this can be integrated to yield 

6(0 = 6(0)6-'^'“* + [ dt' eP'^^'gkit'). (All) 

Assuming that we begin with zero noise, 6(0) = 0; 
Eq. (lAllI) implies that 

(6(06fc(0) = exp -(pfc -P p_fe)t X 

dt'dt" exp {pkf + p-kt''){gk{t')g-k{t"))■ 



dt 


(A12) 
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Using Eq. (|T^ and letting t —>• oo, to obtain the result 
in the stationary state, one finds that if {pk + P-k) > 0, 

hm {ik{t)i-k{t)) = (A13) 

[pk + P-k ) 

Appendix B: Numerical simulation of the 
individual-based model 

The numerical simulations of the individual-based 
model are performed using the Gillespie algorithm [?7l| . 
Before providing the details of our algorithm, we recall 
some basic elements of the process. The state of the 
system at time t is described by N{t) real numbers, 
Xi,i = 1,..., N(t), corresponding to the phenotypes of 
the N{t) individuals present. The probability that indi¬ 
vidual i initiates a reproduction event, given that a birth 
takes place, is 1/N{t) (since we consider no Allee effect, 
the reproduction probability is independent of the den¬ 
sity of suitable mates); the probability that individual i 
chooses individual j to mate, given that individual i is 
reproducing, is proportional to m(xi — Xj) (see Sec. HB: 
finally, the probability density that the newborn individ¬ 
ual has phenotype x, given that individuals i and j are 
reproducing, is given by r{x — {xi -I- Xj)l2) (again, see 
Sec. HI)) . 

There are two possible events (assuming N{t) > 1): 

(i) Death of an individual, with a rate (probability 
per unit of time) equal to ri = / ^{x, 4>)dx = 

9{xi - Xj)/Kf{xi) (see Eqs. ((T|) and dS])). 

(ii) Birth of a new individual, with a rate r 2 = 

/ I3(x, 4>)dx = - Xj)/ - xi) 

[=N{t)) (see Eq.’®). 

Therefore the probability that individual with phenotype 
Xi (that we will denote as individual i) dies, given that a 
death event takes place, is proportional to g{xi — 

Xj)/Kf{x,). 


The numerical simulations are, then, based on the fol¬ 
lowing algorithm: 

1. Set the initial state of the system, that is, the ini¬ 
tial number of individuals, N(to), their correspond¬ 
ing phenotypes, Xi,i = 1,..., N{t), and the initial 
time, t = to- 

2. Compute ri and r 2 , as described earlier. Compute 
the time increment after which the next event takes 
place (At), which is an exponential random variable 
with average (ri -I- 7 - 2 )“^, so it can be computed as 
At = — ln(u)/(ri + r 2 ), with u a pseudo-random 
number with a uniform distribution in the interval 
(0,1). Update the time t = t -|- At. 

3. Establish what type of event takes place. With 
probability ri/(ri -I- r 2 ) a death event takes place; 
go to 4a. Otherwise a birth event takes place; go 
to 4b. 

4a. Establish which individual dies. Choose individ¬ 
ual i at random, with probability proportional to 
'Y/!j=i9{xi — Xj)/Kf{xi). Eliminate individual i. 
Update N, N = N — 1. IfiV = 0 the population 
becomes extinct and the simulation ends. 

4b. Set the phenotype of the new individual. Choose 
individual i uniformly at random to initiate re¬ 
production. Choose individual j to mate, at ran¬ 
dom with probability proportional to m{xi — xj). 
Set the phenotype of the new individual as a: = 
{xi + Xj)/2 + (^, where C is a random variable with 
probability density function given by r(x). Update 
N,N = N + 1 

5. Go to 2 or finish. 

When considering periodic boundary conditions, this has 
to be taken into account when computing the competi¬ 
tion and mating functions as well as the phenotype of the 
new individual. 
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